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Detection  and  classification  of  unexploded  ordnance  based  on  electromagnetic  induction  have  made  tremen¬ 
dous  progress  over  the  last  few  years,  to  the  point  that  not  only  more  realistic  terrains  are  being  considered 
but  also  more  realistic  questions  -  such  as  when  to  stop  digging  -  are  being  posed.  Answering  such  questions 
would  be  easier  if  it  were  somehow  possible  to  see  under  the  surface.  In  this  work  we  propose  a  method  that, 
within  the  limitations  on  resolution  imposed  in  the  available  range  of  frequencies,  generates  subsurface  im¬ 
ages  from  which  the  positions,  relative  strengths,  and  number  of  targets  can  be  read  off  at  a  glance.  The  meth¬ 
od  seeds  the  subsurface  with  multiple  dipoles  at  known  locations  that  contribute  collectively  but 
independently  to  the  measured  magnetic  field.  The  polarizabilities  of  the  dipoles  are  simultaneously  updated 
in  a  process  that  seeks  to  minimize  the  mismatch  between  computed  and  measured  fields  over  a  grid.  In 
order  to  force  the  polarizabilities  to  be  positive  we  use  their  square  roots  as  optimization  variables,  which 
makes  the  problem  nonlinear.  The  iterative  update  process  guided  by  a  Jacobian  matrix  discards  or  selects  di¬ 
poles  based  on  their  influence  on  the  measured  field.  Preliminary  investigations  indicate  a  fast  convergence 
rate  and  the  ability  of  the  algorithm  to  locate  multiple  targets  based  on  data  from  various  state-of-the-art 
electromagnetic  induction  sensors. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Sensing  in  the  electromagnetic  induction  (EMI)  regime  of  frequen¬ 
cies,  which  ranges  from  about  1 0  Hz  to  a  few  hundred  kHz,  has  emerged 
as  an  alternative  to  magnetometry  (Billings,  2004)  and  ground- 
penetrating  radar  (Chen  et  al.,  2001)  in  the  search  for  unexploded 
ordnance  -  the  most  pressing  environmental  problem  faced  by  the 
military  worldwide  (Bilukha  et  al.,  2007;  Sethi  and  Krug,  2000),  and 
one  with  enormous  humanitarian  implications  -  and  in  general  to 
look  for  objects  buried  in  the  ground  (Baum,  1999;  Byrnes,  2009).  The 
process  in  general  consists  of  three  stages:  (a)  Detection  —  inspect  the 
ground  and  take  a  closer  look  at  spots  where  anomalies  are  present; 
(b)  Inversion  —  digest  close-interrogation  data  with  a  forward  model 
that  estimates  the  location  and  orientation  of  the  target  or  targets  that 
cause  the  anomaly  and  synthesizes  an  electromagnetic  signature  for 
each;  and  (c)  Classification  —  distill  electromagnetic  signatures  into 
intrinsic  features  that  allow  the  unambiguous  characterization  of 
targets,  identifying  them  first  and  foremost  as  dangerous  ordnance  or 
innocuous  scrap  and  subsequently  sorting  the  unexploded  munitions 
by  caliber  and  type.  Headway  has  been  made  on  all  three  fronts: 
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concurrent  with  the  design  and  assembly  of  novel  sensors  that  provide 
data  of  unprecedented  quality  and  diversity  has  been  the  development 
of  several  powerful,  fast,  and  accurate  data-processing,  modeling,  and 
classification  algorithms. 

On  the  sensor  side,  the  last  few  years  have  seen  the  appearance  of  in¬ 
struments  that  go  well  beyond  the  upright-dipole-transmitter/vertical- 
component-measuring  architecture  standard  for  beachgoing  metal 
detectors.  Among  cart-based  devices,  the  Time-Domain  Electromag¬ 
netic  Induction  Towed  Array  Detection  System  (TEMTADS)  developed 
by  the  Naval  Research  Laboratory  (Steinhurst  et  al.,  2010)  replicates 
this  coplanar/concentric  transmitter/receiver  paradigm  on  a  5  x  5  grid, 
producing  625  time  histories  per  instrument  location,  while  the 
Geometries  MetalMapper  (Prouty,  2009a;  Prouty,  2009b)  uses  fewer 
modules  but  arranges  both  transmitters  and  receivers  to  explore  the 
full  three-dimensional  vector  profile  of  the  field.  Hand-based  instru¬ 
ments  like  the  Man  Portable  Vector  (MPV)  sensor  (Barrowes  et  al., 
201 1 ;  Fernandez  et  al.,  201 1 )  have  also  incorporated  multiple  three- 
axis  vector  receivers  as  part  of  the  design  and  allow  sensing  in  uneven 
or  wooded  terrain.  These  next-generation  instruments  are  ultrawide- 
band  (i.e.,  they  nearly  span  the  complete  range  of  EMI  frequencies) 
and  work  either  in  the  frequency  domain  (Fernandez  et  al.,  2009; 
Won  et  al.,  1997;  Won  et  al.,  2001)  or  in  the  time  domain,  where  they 
can  take  readings  lasting  up  to  25  ms  after  transmitter  shutoff. 
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On  the  analysis  side,  several  new  procedures  have  either  honed  or 
extended  the  standard  model  that  assumes  every  target  to  be  a  single 
point  dipole  (Baum,  1999;  Pasion  and  Oldenburg,  2001).  The  EMI- 
response  problem  has  been  solved  analytically  for  spheroids  (Ao  et 
al.,  2002;  Barrowes  et  al.,  2004;  Braunisch  et  al,  2001;  Chen  et  al, 
2007)  and,  partially,  for  ellipsoids  (Grzegorczyk  et  al.,  2008).  Some 
semi-analytic  models  expand  the  fields  in  geometrically  meaningful 
modes  and  solve  the  problem  mode  by  mode  (Fernandez  et  al., 
2008;  Shubitidze  et  al.,  2004;  Sun  et  al.,  2005).  Other  methods  replace 
the  point  dipole  by  surface  distributions  of  magnetic  charge 
(Shubitidze  et  al.,  2007)  or  dipole  moment  (Shubitidze  et  al.,  2010a) 
spread  over  a  surface  that  encloses  the  target;  more  recent  versions 
use  volumetric  dipole  distributions  (Shubitidze  et  al.,  2010b).  Other 
possibilities  include  direct  optimizations  based  on  analytical  methods 
(Zhang  and  Liu,  2001),  model-based  procedures  (Miller  et  al.,  2001), 
and  statistical  approaches  (Tantum  and  Collins,  2001). 

The  resulting  combination  of  hardware  and  software  has  been 
subjected  to  a  battery  of  tests  of  increasing  realism  in  places  like  the 
former  Camp  Sibert  in  Alabama  (Fernandez  et  al.,  2010b;  Shamatava 
et  al.,  2009;  Shubitidze  et  al.,  2010a),  the  Aberdeen  Proving  Ground 
(APG)  in  Maryland  (Fernandez  et  al.,  2010a;  Shubitidze  et  al.,  2009), 
the  former  Camp  San  Luis  Obispo  (SLO)  in  California  (Prouty,  2009a; 
Prouty,  2009b;  Shamatava  et  al.,  2010),  and  Camp  Butner  in  North 
Carolina  (Shamatava  et  al.,  2011).  The  results  have  been  consistently 
encouraging,  and  for  that  reason  have  opened,  rather  than  closed,  the 
need  for  more  sophisticated  methods:  the  UXO  community  must  now 
consider  ever  more  difficult  questions  based  on  acute  real-world 
problems.  For  example,  the  need  to  decontaminate  littoral  terrain 
(SanFilipo  et  al.,  2005)  has  prompted  the  study  of  the  effects  of  a 
conductive  embedding  medium  on  EMI  sensing.  The  high  densities  of 
both  clutter  and  legitimate  targets  present  in  many  UXO  sites  has  led 
to  the  problem  of  inverting  for  the  locations  and  electromagnetic  signa¬ 
tures  of  multiple  targets  present  in  the  spatial  range  of  a  sensor 
(Grzegorczyk  et  al.,  2009;  Song  et  al.,  2009;  Grzegorczyk  et  al.,  2011) 
in  order  to  approximate  either  multi-UXO  configurations  or  UXO  in 
the  presence  of  a  few  bulky  clutter  items.  Indeed,  it  is  still  not  clear 
how  to  deal  consistently  with  clutter:  questions  like  how  large  a  piece 
of  clutter  should  be  to  be  inverted  as  an  additional  target,  or  what 
density  and  size  distribution  a  cluster  of  clutter  items  should  have  to 
be  treated  as  noise  rather  than  a  collection  of  individual  objects,  are 
not  at  all  settled. 

EMI  frequencies  (again,  between  tens  of  Hz  and  hundreds  of  kHz) 
are  useful  for  subsurface  sensing  because  the  ground,  with  typical  con¬ 
ductivity  a-  0.1-1000  mS/m  (though  the  high  end  of  the  range  is  very 
rare),  permittivity  £~ 50-177  pF/m,  and  permeability  /i^/io-10-6 
H/m,  has  a  skin  depth  of  the  order  of  tens  of  meters  at  least  and,  there¬ 
fore,  is  transparent  in  the  range  of  interest  (Keller,  1987;  Milsom  and 
Eriksen,  2011).  (The  values  just  quoted  do  not  exhaust  the  very  ample 
ranges  in  electromagnetic  properties  exhibited  by  soils  (Butler,  2003), 
for  example  in  places  like  Hawai‘i,  but  are  representative  of  the  majority 
of  non-littoral  soils  of  interest  to  the  UXO  community.)  This  frequency 
range,  however,  is  characterized  by  a  very  low  resolution,  which, 
together  with  the  fact  that  sensing  must  be  done  from  a  distance, 
above  ground,  means  that  one  cannot  form  any  reasonably  clear  images 
of  the  subsurface  objects.  On  the  other  hand,  the  available  resolution 
does  suffice  to  locate  objects  in  three  dimensions  and  to  resolve  targets 
separated  by  a  few  centimeters.  This  ability  -  which  depends  on  the 
sensor  employed  and  on  the  sizes,  orientations,  and  composition  of 
the  relevant  targets  -  opens  the  possibility  of  generating  subsurface 
images  to  aid  in  the  UXO  remediation  process.  In  this  paper  we  explore 
this  possibility  and  propose  a  methodology  to  carry  it  out. 

In  current  practice,  responsibility  for  estimating  the  number  of 
targets  producing  a  given  anomaly  falls  upon  the  inversion  stage, 
which  typically  means  that  a  time-consuming  and  computationally 
expensive  inversion  procedure  has  to  be  run  either  assuming  a 
sequential  number  of  targets  (one,  two,  three,  etc.),  or  assuming  an 


upper  bound  on  the  number  of  targets,  in  which  case  the  redundant 
dipoles  either  cluster  around  expected  positions  or  take  very  small 
amplitudes  that  immediately  expose  them  as  irrelevant.  An  increasing 
number  of  targets  increase  the  number  of  model  parameters  to  be 
inverted  and  riddles  the  search  space  with  local  minima,  making  the 
inversion  less  reliable  and  costlier  in  CPU  time  —  a  serious  shortcoming 
in  a  field  where  real-time  analysis  is  desirable. 

Having  a  subsurface  image,  on  the  other  hand,  would  immediately 
give  operators  a  set  of  visual  cues  to  guide  their  intuition.  It  would 
provide  an  initial  estimate  of  the  number  of  targets,  with  which  it 
would  be  possible  to  run  the  inversion  routine  just  once;  the  image 
would  also  yield  positional  information  that  could  accelerate  the  inver¬ 
sion  by  providing  a  sound  initial  guess  for  the  optimization.  The  images, 
moreover,  could  allow  essentially  instantaneous  anomaly  prioritization 
by  allowing  the  user  to  pick  out  the  most  prominent  anomalies,  and 
could  also  help  distinguish  between  large,  deep  targets  and  small,  shal¬ 
low  ones.  Furthermore,  an  “imaging”  survey  at  the  beginning  of  the 
remediation  would  let  operators  make  educated  guesses  about  the 
site  -  i.e.,  the  estimated  number  and  density  of  dangerous  targets, 
the  clutter  distribution,  the  regions  to  which  more  attention  must  be 
devoted  -  and  thus  make  more  efficient  use  of  their  time  and  resources. 
In  all,  an  image,  as  opposed  to  a  beep  or  a  number,  makes  the  detection 
of  UXO  more  visual  and  user-friendly,  and  eventually  more  reliable 
because  it  enables  the  user  to  cooperate  with  the  apparatus. 

The  work  presented  in  (Miller  et  al.,  2010)  reports  promising  results 
in  this  direction  found  using  harmony  search,  a  method  in  which  “trials” 
with  a  preset  number  of  sources  are  sequentially  generated  -  either 
afresh  or  by  making  slight  changes  to  existing  configurations  -  and 
accepted  or  discarded  depending  on  their  “success”  —  i.e.,  the  closeness 
of  their  predictions  to  measured  values.  The  authors  manage  to  properly 
resolve  multi-target  data  taken  with  TEMTADS  at  the  Blossom  Point 
Army  Research  Facility  in  Maryland  and  at  APG.  The  method  we 
propose  here  is  fundamentally  different  in  that  it  is  not  stochastic.  As 
in  (Miller  et  al.,  2010),  we  (a)  seed  the  region  below  the  surface  with 
a  distribution  of  point  dipoles,  each  of  which  contributes  to  the 
measured  magnetic  field,  and  (b)  iteratively  update  the  dipoles' 
polarizabilities  -  that  is,  their  dipole  moments  but  dividing  out  the  pri¬ 
mary  field  that  induces  them-so  as  to  minimize  the  mismatch  between 
computed  and  measured  fields  over  a  grid,  but  here  the  similarities  end. 
Instead  of  looking  for  the  minimum  sufficient  number  of  sources,  we 
look  for  the  sources  that  have  the  most  influence  on  the  measured 
data  and  iteratively  “zoom  in”  on  them.  Instead  of  avoiding  local 
minima  by  occasionally  rousting  the  solution  with  a  random  ripple, 
we  use  a  Jacobian  to  guide  the  search  downhill.  We  also  force  the 
positivity  of  the  polarizabilities  by  using  their  square  roots  as  the  pa¬ 
rameters  to  be  estimated;  this  imposition  makes  the  problem  nonlinear, 
but  rids  the  results  of  trivial  solutions  and  improves  convergence  to  the 
global  minimum. 

2.  Procedure 

Our  approach  is  based  on  the  discretization  of  the  subsurface  into 
dipoles  that  contribute  collectively  and  independently  to  the  measured 
magnetic  field.  The  locations  of  the  dipoles  are  known  by  construction, 
but  not  their  polarizabilities,  and  the  task  of  the  algorithm  is  to  compute 
the  polarizabilities  that  minimize  the  difference  between  measured  and 
computed  fields.  The  polarizability  amplitudes  that  result  should  be 
high  near  the  true  locations  of  the  targets  and  low  everywhere  else, 
allowing  the  generation  of  images.  There  are  three  polarizability 
elements  per  dipole,  each  of  which  can  either  be  used  by  itself  to  synthe¬ 
size  an  image  -  yielding  not  just  estimates  of  a  target's  location  but  also 
of  its  orientation  -  or  combined  with  the  others  to  form  a  composite. 

As  an  initial  proof  of  concept  we  study  multitarget  ensembles 
arranged  in  such  a  way  that  it  is  possible  to  find  a  vertical  plane  that 
simultaneously  cuts  either  through  or  close  to  all  of  the  targets.  We 
then  take  this  plane  to  define  a  common  and  known  y-coordinate, 
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distribute  a  two-dimensional  grid  of  Nd  dipoles  on  it,  and  make  it  our 
task  to  synthesize  images  that  let  us  determine  the  unknown  locations 
of  the  targets  starting  from  field  data.  The  restriction  to  a  vertical  2D 
domain  of  sources  and  the  placement  of  the  sources  on  a  regular  grid 
have  been  adopted  in  this  paper  for  ease  of  implementation  and 
speed  of  execution,  but  neither  is  a  limitation  of  the  method,  which 
we  expect  to  be  directly  generalizable  to  3D  configurations.  The  issues 
that  naturally  arise  in  such  generalizations  -  computational  burden, 
speed,  and  memory  requirements  -  appear  as  potential  challenges  that 
need  to  be  addressed;  this  task  we  postpone  to  future  considerations. 

The  source  dipoles,  though  constrained  to  lie  in  a  plane,  respond  in 
a  fully  three-dimensional  manner  to  the  primary  field  of  the  sensor 
through  the  standard  expression  (Van  Bladel,  1964) 


Hd(r„) 


1/3 r(qm)__ 
4rr|r |3  V  |r|2 


(1) 


where  r  points  from  the  dipole  to  the  observation  point  r0  and  the  di¬ 
pole  moment 

m  =  xmx  +  ymy  +  zmz 

=  x/3xHpxr  +ypyri?  +  zf3zHpzr, 

with  Hpr  the  primary  field  at  the  location  of  the  dipole.  The  polariz¬ 
abilities  £e{x,j/,z},  are  in  general  elements  of  a  symmetric  3x3 
tensor.  We  force  the  principal  axes  of  all  the  dipoles  to  coincide 
with  those  of  the  sensor,  resulting  in  a  diagonal  polarizability  tensor 
with  only  three  nonzero  values.  We  thus  set  out  to  solve 


We  use  a  Gauss-Newton  method  (Press  et  al.,  1992)  for  the  task,  regu¬ 
larized  with  a  Tikhonov  term  in  order  to  give  preferential  weight  to  so¬ 
lutions  with  smaller  norms  (Aster  et  al.,  2005).  The  3Ndxl  solution 
vector  x,  equivalent  to  the  elementwise  square  root  of  y,  is  formally 
the  limit  as  i  — ► 00  of  an  iteration-dependent  x1.  The  initial  value  x°  is  cho¬ 
sen  at  random,  though  making  reasonable  assumptions  (e.g.,  that  the  po¬ 
larizabilities  have  roughly  the  same  order  of  magnitude),  and  each 
subsequent  update  to  the  vector  comes  from  solving  the  normal  equation 

xi+1  =  x'  +  (]  T  J  +  Al)  _1  J  T  AH  (x),  (7) 


where  the  /-dependent  N  •  Ncmpx3Nd  Jacobian 


j  = 


dH^/dx1' 

dfi2/dxl 


_dHN/dxl. 


(8) 


guides  the  search  by  approximating  the  Hessian  (Gill  et  al.,  1981 ).  The  de¬ 
rivatives  with  respect  to  the  square  roots  of  the  polarizabilities  are 
straightforward  to  compute,  given  the  simple  dependence  of  Eq.  (1)  on 
the  {/3}.  The  Tikhonov  regularization  parameter  \  (Hansen,  2010)  is  ini¬ 
tially  chosen  so  that  the  condition  number  of  the  regularized  approximate 

Hessian  matrix  (jTJ  +  A/)  is  no  less  than  10“ 10,  and  then  adjusted  at 
every  iteration  using 


S  t.  {/\}  >0  v$, 


where  we  have  enclosed  the  /3^  in  curly  brackets  to  emphasize  that 
the  minimization  involves  all  Nd  seeded  dipoles  and  have  strung  to¬ 
gether  the  measured  fields  in  the  vector 


j  jHlcdS _ 


H? 


(4) 


of  length  N-  NcmpE  (Nioc  •  NTx  •  NrJ  •  Ncmp,  where  Nioc  is  the  number 
of  sensor  locations  at  which  data  are  collected,  NTx  and  Nrx  are  re¬ 
spectively  the  number  of  transmitters  and  receivers  on  the  sensor, 
and  Ncmp  is  the  number  of  measured  field  components.  We  also  have 
made  explicit  the  fact  that  the  polarizabilities  are  always  positive 
(Pasion  et  al.,  2008);  we  impose  this  constraint  in  the  minimization  by 
using  the  square  roots  of  the  polarizabilities,  =  V  as  inversion  pa¬ 
rameters: 


min 

It's} 


+  M2 


=  min  (| |AH (y)||2  +  ||Ay||2) , 

M 

where 


(5) 


y /  = 

,  y  = 

"yi  ‘ 

Wz,)\ 

Vn,. 

the  Green  dyads  G  incorporate  the  dipole  contributions  from  Eq.  (3) 
with  the  primary  field  factored  in,  and  AH  stands  for  the  mismatch  be¬ 
tween  modeled  and  measured  signals.  Using  7 $  as  the  variables  instead 
of  j8^  makes  the  problem  nonlinear,  one  that  has  to  be  solved  iteratively. 


where  £  is  adjustable  (we  use  g=  10“ 4  throughout  this  paper).  At  the 
outset  the  matrices  are  ill-conditioned,  and  the  regularization  term 
helps  stabilize  the  inversion;  as  the  iteration  progresses,  the  field  mis¬ 
match  AH  tends  to  zero  and  A  becomes  ever  less  important  compared 
to  the  data. 

Once  the  process  converges  we  look  for  the  region  with  the  highest 
concentration  of  polarizability  amplitude  and  zoom  in  on  it  with  the 
aim  of  generating  more  finely  detailed  images.  Restricting  our  attention 
to  the  smaller  computational  space,  we  distribute  a  similar  grid  of 
dipoles  and  perform  a  new  optimization.  This  refocusing  helps  us 
keep  a  moderate  number  of  source  dipoles  and  makes  for  a  more  effi¬ 
cient  algorithm;  let  us  also  emphasize  that  no  additional  sensor  data 
are  required.  The  process  is  repeated  until  we  obtain  a  “polarizability 
map,”  with  high  concentrations  of  polarizability  being  indicative  of 
target  location.  At  this  stage,  the  refocusing  is  performed  manually,  by 
eye,  but  it  is  not  hard  to  envision  how  the  process  could  be  automated. 

3.  Results 

We  have  performed  initial  validating  studies  of  our  algorithm  using 
both  synthetic  and  measured  data.  The  synthetic  data  incorporate  the 
geometries  of  current  sensors  and  targets  with  realistic  signatures  and 
locations. 

3.1.  Synthetic  scenarios 

Initially  we  compare  the  images  that  can  be  obtained  when  the  algo¬ 
rithm  has  access  to  full  three-dimensional  magnetic-field  data  and 
when  only  the  vertical  component  is  available.  For  definiteness  we 
consider  a  TEMTADS-like  configuration:  we  generate  images  starting 
from  a  single  data  shot  as  collected  by  25  receivers  distributed  on  a 
5x5  rectangular  uniform  grid,  with  each  receiver  separated  from  its 
neighbors  center-to-center  by  40  cm.  The  sensor  is  assumed  to  lie 
centered  at  the  origin  and  flat  on  the  ground.  To  simulate  unfavorable 
measurement  conditions,  we  take  only  a  subset  of  the  data  that  can  pos¬ 
sibly  be  collected  by  the  instrument:  though  each  TEMTADS  receiver  is 
concentric  with  a  transmitter,  we  use  the  data  collected  when  only  the 
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center  transmitter  is  fired.  The  targets  are  two  point  dipoles  placed  at 
x  =  ±  10  cm  at  a  depth  of  50  cm  below  the  sensor  center  and  oriented 
along z  with  unit  magnitude  (hence  /3X  =  /3y  =  0,  /3Z  =  1  for  each  dipole). 
A  set  of  50x50  source  dipoles  are  distributed  on  a  two-dimensional 
vertical  grid  on  the  y  =  0  plane  spanning  2  m  (the  side  length  of  TEM- 
TADS)  along  the  x-direction  and  going  from  20  cm  to  1  m  in  depth. 
Figs.  1  and  2  respectively  report  the  results  found  using  all  three  compo¬ 
nents  of  H  and  using  Hz  only.  The  figures  present  the  normalized  polar¬ 
izability  magnitude  of  the  jth  dipole,  computed  using 

P,=Pj/nux(Pj),  0 j  =  y/&j+%j+gj,  j=\,..,Nd.  (10) 

The  bottom-right  panels  in  the  figures  illustrate  the  typical  conver¬ 
gence  rate  of  our  algorithm  by  showing  how  the  normalized  error 
evolves  as  a  function  of  the  iteration  number.  Note  that  in  Figs.  1 
and  2  the  pixel  amplitudes  of  images  corresponding  to  {/3X}  and  {j8y} 
are  negligible  compared  to  those  corresponding  to  {jBJ,  as  expected, 
and  are  not  shown.  The  colors  in  the  figures  range  from  blue  to  red,  cor¬ 
responding  to  respective  amplitudes  of  zero  and  one,  with  a  similar 
format  adopted  for  all  subsequent  figures.  The  figures  illustrate  the 
fact  that  polarizability  amplitudes  are  much  stronger  close  to  the  target 
locations  than  everywhere  else,  where  they  are  essentially  zero.  Physi¬ 
cal  information  on  dipole  strength  can  be  obtained  by  integrating  the 
images. 

Though  in  both  cases  the  targets  are  properly  resolved  and  found 
at  locations  not  far  from  the  ground  truth  -  denoted  in  the  figures  by 
white  circles  -  it  is  not  surprising  that  the  objects  are  easier  to  discern 
in  the  former  situation:  just  one  zoom  is  necessary  to  resolve  properly 


the  targets  in  Fig.  1 ,  while  they  remain  poorly  separated  in  the  case  of 
Fig.  2. 

We  remind  the  reader  that  this  hypothetical  exercise  aims  to  study 
the  effect  of  using  one  vs.  three  components  of  the  magnetic  field  and 
is  not  expected  to  be  generalized  to  the  actual  TEMTADS  sensor 
array,  which  provides  625  data  points  rather  than  the  25  we  use  here. 
A  realistic  TEMTADS  case  is  considered  below. 

As  a  second  example,  we  consider  a  situation  in  which  a  weaker 
dipole  (which  could  represent  a  clutter  item)  is  buried  at  a  shallow 
depth  above  a  deeper  yet  stronger  target  that  could  stand  for  UXO. 
This  is  a  common  occurrence  in  the  field,  and  one  that  time  and 
again  has  proved  difficult  to  analyze.  We  use  the  same  TEMTADS- 
like  configuration  from  the  previous  example  and  use  all  components 
of  the  magnetic  field.  The  targets  are  placed  10  cm  to  the  right  of  the 
simulated  sensor  in  the  x-direction  and  centered  in  y;  the  “clutter” 
target  is  at  a  depth  of  20  cm;  the  “UXO”  target  is  30  cm  deeper  and 
its  dipole  moment  has  twice  the  magnitude.  Both  objects  are  oriented 
alongz.  The  source  dipoles  are  on  the  same  initial  grid  as  in  the  previous 
example.  The  results,  which  appear  in  Fig.  3,  again  reveal  that  the  algo¬ 
rithm  can  identify  the  region  in  which  targets  are  present  and  can 
distinguish  the  shallow  object  from  the  deeper,  stronger  one. 

3.2.  Experimental  configurations 

We  subsequently  illustrate  the  performance  of  our  algorithm  on 
actual  measurements  taken  with  three  different  state-of-the-art  EMI 
sensors:  the  TEMTADS  array  (Steinhurst  et  al.,  2010),  the  MPV  sensor 
(Fernandez  et  al.,  2011;  George  and  George,  2007),  and  the  Metal- 
Mapper  system  (Prouty,  2009a).  Each  of  these  sensors  presents  a 
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Fig.  1.  Image  based  on  synthetic  data  generated  over  a  5  x  5  grid  similar  to  that  of  TEMTADS  but  using  only  one  transmit  coil.  Two  dipoles  of  similar  moment  are  placed  at  the  locations  in¬ 
dicated  by  the  white  circles.  All  vector  components  of  the  polarizabilities  are  added  pixel-to-pixel  to  generate  the  image.  Pixel  amplitudes  are  normalized  between  0  (blue)  and  1  (red).  Top- 
left:  First  stage  with  large  initial  region.  Top-right:  Second  stage  and  refocusing  within  the  red  rectangle  of  (A).  Bottom-left:  Third  stage  and  refocusing  within  the  red  rectangle  of  (B).  Bottom- 
right:  Normalized  error  as  a  function  of  the  iteration  number  for  (A). 
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Fig.  2.  Same  as  Fig.  1  but  using  only  the  z-component  of  the  magnetic  field.  It  is  seen  that  two  refocusing  stages  are  needed  to  resolve  the  two  sources.  Pixel  amplitudes  are  nor¬ 
malized  between  0  (blue)  and  1  (red).  Top-left :  First  stage  with  large  initial  region.  Top-right :  Second  stage  and  refocusing  within  the  red  rectangle  of  (A).  Bottom-left:  Third  stage 
and  refocusing  within  the  red  rectangle  of  (B).  Bottom-right:  Normalized  error  as  a  function  of  the  iteration  number  for  (A). 


different  transmit/receive  configuration  that  is  unique  enough  to 
warrant  independent  processing  with  our  algorithm. 

3.2.1.  TEMTADS  data 

The  TEMTADS  array,  already  described  above,  is  composed  of  a  5  x  5 
square  grid  of  parallel  coplanar  flat  transmitter/receiver  pairs.  A  total  of 
625  measurements,  corresponding  to  each  transmitter/receiver  combi¬ 
nation,  are  available  for  analysis  at  each  sensor  location,  and  the  40-cm 
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Fig.3.  Normalized  polarizability  map  obtained  from  synthetic  data  generated  over  a 
TEMTADS-like  5x5  grid  with  a  single  transmit  coil.  Two  dipoles  are  placed  on  top  of 
each  other  at  (10,0,  —  20)  cm  and  (10,0,  —  50)  cm,  as  indicated  by  the  white  circles.  The 
shallow  dipole  has  half  the  strength  of  the  deeper  one.  The  images  are  presented  as  a  sin¬ 
gle  polarizability  magnitude  obtained  from  Eq.  (10). 


neighbor-to-neighbor  separation  generates  data  of  enough  diversity  in 
viewpoint  to  make  up  for  the  fact  that  only  the  z-component  of  the 
field  is  measured. 

Our  TEMTADS  example  has  three  targets  on  the  x-z  plane:  a  mortar 
buried  60  cm  below  the  origin  and  pointing  along  x,  a  baseplate  at  a 
depth  of  44  cm  and  50  cm  to  the  left  of  the  mortar,  and  a  nosepiece 
buried  29  cm  below  ground  and  20  cm  to  the  left  of  the  baseplate  (see 
Fig.  4).  In  a  previous  paper  (Grzegorczyk  et  al.,  201 1 )  we  analyzed  this 
same  configuration,  estimating  the  locations,  orientations,  and  entire 
time-dependent  polarizabilities  of  all  three  targets  using  a  full  inver¬ 
sion/discrimination  algorithm  but  knowing  a  priori  the  number  of 
buried  targets.  In  the  present  work  we  make  no  assumptions  about 
the  number  of  targets  but  attempt  to  estimate  it  from  the  image  itself. 
In  addition,  we  focus  on  target  location  and  not  on  time-dependent 
signatures,  and  for  that  reason  run  the  imaging  algorithm  on  a  single 
time  channel,  typically  an  early  one  with  high  signal-to-noise  ratio. 

The  subsurface  in  this  example  is  seeded  with  40  x  40  dipoles  over 
a  uniformly  distributed  rectangular  grid  on  the  y  =  0  plane  spanning 
between  —  1  m  and  +1  m  in  the  x-direction  and  between  10  cm 
and  1  m  in  depth.  The  uniform  dipole  distribution  is  not  a  require¬ 
ment  of  the  algorithm  but  merely  reflects  the  fact  that  we  assume 
no  previous  knowledge  about  the  configuration  beyond  the  shared 
y-location.  The  algorithm  is  run  with  10  iterations  and  produces  the 
results  shown  in  Fig.  4.  The  superposed  white  circles  indicate  the 
ground  truth.  Overall  we  see  a  good  match  between  the  actual  target 
locations  and  the  peaks  in  polarizability;  the  background  is  seen  to  be 
void  of  spurious  sources.  We  note  that  the  nosepiece,  despite  being 
the  target  closest  to  TEMTADS,  is  the  most  difficult  to  find;  we  attribute 
the  inaccuracy  to  the  nosepiece's  being  almost  at  the  edge  of  the  sensor 
and  therefore  interrogated  from  relatively  few  angles. 
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Fig.  4.  Configuration  and  resulting  normalized  polarizability  map  obtained  from  a 
three-target  TEMTADS  measurement.  The  targets  are  as  follows:  mortar  at  (0,0,— 60)  cm, 
baseplate  at  (  —  50,0,-44)  cm,  and  nose  piece  at  (-70,0,-29)  cm.  Axes  dimensions  of 
the  image  are  in  meters. 


3.2.2.  MPV  data 

As  its  name  indicates,  the  MPV  (described  in  more  detail  else¬ 
where  (Fernandez  et  al.,  2011))  is  a  man-portable  and  therefore 
much  smaller  sensor.  It  has  a  single  horizontal  transmit  direction  for 
the  primary  field  but  captures  all  three  components  of  the  secondary 
field  at  five  locations,  one  of  them  at  a  different  height  than  the 
others.  The  field  of  view  is  consequently  much  smaller  than  that  of 
TEMTADS,  though  this  limitation  (as  well  as  the  lack  of  primary- 
field  diversity)  could  be  overcome  by  waving  the  instrument  around 
and  tilting  it  above  an  anomaly.  This,  however,  is  not  attempted  here. 

Our  MPV  example  has  two  40-mm  UXO  within  the  area  scanned 
by  the  sensor,  located  at  the  same  depth  of  40  cm  and  separated  by 
a  lateral  distance  also  of  40  cm.  The  data  were  collected  by  Sky  Re¬ 
search  personnel  in  Hanover,  New  Hampshire,  using  a  5  x  5  uniform 
square  grid  of  side  80  cm  and  a  point-to-point  separation  of  20  cm. 
The  algorithm  is  run  using  a  uniform  40  x  40  dipole  distribution  seed¬ 
ed  on  the  x-z  plane  and  confined  to  a  lateral  region  between  —  50  cm 
and  +  50  cm  in  x  and  between  10  cm  and  50  cm  in  depth.  The  result, 
shown  in  Fig.  5,  reveals  the  presence  of  two  targets  at  locations  that 
match  the  ground  truth  (white  circles)  well. 

3.2.3.  MetalMapper  data 

As  a  final  example,  we  consider  data  collected  by  the  MetalMapper 
in  a  test  mode  at  Camp  San  Luis  Obispo  (Prouty,  2009a;  Prouty, 
2009b).  The  instrument  has  three  mutually  perpendicular  transmitter 
coils  and  seven  receiver  cubes  that  measure  all  three  vector  compo¬ 
nents  of  the  secondary  field.  The  dataset  is  therefore  rich  in  information, 
though  geometrically  more  restricted  than  those  achievable  with  TEM¬ 
TADS.  The  targets  are  two  copper  rings,  each  of  19  cm  and  21  cm  inner 
and  outer  diameters  and  otherwise  thin.  One  is  buried  flat  under  the 
first  quadrant  of  the  sensor;  the  other  is  placed  vertically  under  the 
fourth  quadrant  (see  Fig.  6).  The  geometry  of  the  latter  target  makes  it 
difficult  to  induce  eddy  currents  in  it  with  a  less  flexible  instrument. 

The  algorithm  employs  40x40  dipoles  uniformly  distributed  be¬ 
tween  —  0.7  m  and  +  0.7  m  in  x  and  between  —  0.6  m  and  —  0.1  m  in 
z.  The  image  of  Fig.  7  results  after  ten  iterations.  As  we  state  at  the  be¬ 
ginning  of  Section  2,  the  procedure  produces  three  images,  one  along 


X  axis  [m] 

Fig.  5.  Normalized  polarizability  map  obtained  from  a  two-target  MPV  measurement. 
The  buried  objects  are  two  identical  40-mm  projectiles  located  at  (0,0, —  40)  cm  and 
(40,0,-40)  cm.  In  this  case  there  is  no  need  to  “refocus”  the  algorithm. 


each  canonical  direction.  The  thin  rings  will  tend  to  be  visible  only 
along  the  direction  where  it  is  possible  to  induce  currents  in  them, 
which  means  that  each  of  them  would  appear  in  only  one  of  the  images. 
Fig.  7  has  been  obtained,  as  elsewhere  in  this  paper,  by  computing  the 
normalized  polarizability  magnitude  of  (10).  Both  targets  are  seen  to 
be  properly  located  when  compared  to  estimates  from  two  indepen¬ 
dent  full  inversions  (which  again  appear  as  white  circles).  We  do  note 
that  the  method  overestimates  the  depth  of  the  deeper  target  by  a 
few  centimeters. 

4.  Conclusion 

In  this  paper  we  have  presented  a  method  that  could  be  used  to 
generate  subsurface  images  and  locate  unexploded  ordnance  using 
already  available  EMI  sensors.  As  validating  examples  we  studied  sever¬ 
al  multi-target  configurations,  some  synthetically  generated  and  others 
measured  by  the  TEMTADS,  MPV,  and  MetalMapper  next-generation 
systems.  The  accuracy  of  our  results  suggests  that  the  method  has  the 
flexibility  required  to  process  varied  datasets  adequately,  including 
those  with  scalar/vector  transmitters  and/or  scalar/vector  receivers. 
Moreover,  our  algorithm  uses  the  data  already  provided  by  existing 
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Fig.  6.  Schematic  representation  of  the  z  MetalMapper  transmitting  coil,  the  seven  vector 
receivers,  and  the  positions  of  the  horizontal  and  vertical  rings.  The  other  transmitters  are 
perpendicular  to  the  x  and  y  directions  and  are  centered  56  cm  above  the  plane  of  the 
sheet. 
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Fig.  7.  Two-target  normalized  polarizability  map  obtained  from  the  MetalMapper  data. 


sensors  and  requires  nothing  more.  We  believe  that  this  method  could 
be  incorporated  into  current  systems  and  in  that  capacity  prove  benefi¬ 
cial  to  the  field  teams  in  charge  of  locating  and  eventually  unearthing 
UXO. 

The  proof  of  concept  presented  in  this  paper  has  concentrated  on  2D 
situations,  with  all  targets  sharing  a  known  y-coordinate.  The  generali¬ 
zation  to  3D  imaging  is  conceptually  straightforward,  given  that  the 
method  is  already  three-dimensional  in  nature,  but  the  expected 
increase  in  computational  needs  poses  challenges  that  may  require 
investigation  of  acceleration  methods. 

Other  limitations  of  the  model  also  need  to  be  investigated  quantita¬ 
tively  before  the  technology  is  mature  enough  for  field  deployment.  In 
particular,  resolvability  with  respect  to  realistic  levels  of  noise  in  EMI 
data  needs  to  be  characterized  in  order  to  quantify  how  small,  how 
deep,  and  (in  multi-target  cases)  how  closely  clustered  subsurface 
objects  can  be  and  still  be  resolvable. 

Results  are  expected  to  vary  not  only  with  the  level  of  noise  asso¬ 
ciated  with  EMI  measurements  and  instrument  positioning,  but  also 
with  the  accuracy  with  which  the  background  can  be  canceled.  In 
our  experience,  real-field  situations  are  very  diverse,  and  resolvability 
should  be  studied  on  a  case-by-case  basis;  for  that  reason  we  decided 
against  adding  synthetic,  one-size-fits-all  noise  to  the  synthetic  cases 
present  in  this  paper.  Finally,  it  would  be  desirable  to  incorporate  into 
the  model  the  time  evolution  of  signals,  given  its  potential  usefulness 
as  a  discrimination  tool.  Such  considerations  represent  the  next  steps 
to  take  based  on  the  results  we  have  presented  here. 
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